#####################################################################################################################################################
#Purpose: To convert plink file to annovar input file
#Date: 11-14-2012
#####################################################################################################################################################

use Getopt::Long;

&Getopt::Long::GetOptions(
'INPUT=s' => \$inputfile,
'SAMTOOLS=s' => \$samtools,
#'chr_end_cord=s' => \$chr_end_cord
'OUTPUT=s' => \$outputfile,
'REFDIR=s' => \$refdir
);
if($inputfile eq "" || $samtools eq "" || $outputfile eq "" || $refdir eq "")
{

	die "missing arguments\n USAGE : perl perl_prep_file_annovar.pl -INPUT <INPUT PLINK BIM FILE> -SAMTOOLS <PATH TO SAMTOOLS BIN> -OUTPUT <OUTPUT FILE> -REFDIR <PATH TO UCSC REF FASTA FILES DIR>\n";
        #die "missing arguments\n USAGE : perl perl_prep_file_annovar.pl -INPUT SNP.bim -SAMTOOLS /projects/bsi/bictools/apps/alignment/samtools/samtools-0.1.18/samtools -OUTPUT SNP1.bim -REFDIR ./ucsc_hg18/\n";
}

chomp($inputfile);
chomp($samtools);
chomp($outputfile);
#perl perl_prep_file.pl SNP.bim /projects/bsi/bictools/apps/alignment/samtools/samtools-0.1.18/samtools SNP1.bim ./ucsc_hg18/
chomp($refdir);
print "Here are the arguments INPUT : $inputfile\n SAMTOOLS : $samtools\n OUTPUT : $outputfile\nREFDIR : $refdir\n";
#die;
open(GNT,$inputfile) or die "no file inputfile found\n";
open(WRGNT,">$outputfile") or die " not able to write file $outputfile found\n";
while(<GNT>)
{
	chomp($_);
	@gnt=split("\t",$_);
	if($gnt[0] ==23)
	{
		$gnt[0] ="X";
	}
	if($gnt[0] ==24)
        {
                $gnt[0] ="Y";
        }
	if($gnt[0] ==26)
        {
                $gnt[0] ="M";
        }
	$sys=`$samtools faidx $refdir/chr$gnt[0].fa chr$gnt[0]:$gnt[3]-$gnt[3]`;
	#@sys=split($sys,"\n";
	$sys =~ s/^.*\n//g;
	$sys =~ s/\n//g;
	$sys_tr =$sys;
	$sys_tr =~ tr/ATGCatgc/TACGtacg/;
	if(uc($gnt[4]) eq uc($sys) || uc($gnt[4]) eq uc($sys_tr)) 
	{
		print WRGNT "$gnt[0] $gnt[3] $gnt[3] $gnt[4] $gnt[5] comments: $gnt[1]\n";
	}
	elsif(uc($gnt[5]) eq uc($sys) || uc($gnt[5]) eq uc($sys_tr))
	{
		print WRGNT "$gnt[0] $gnt[3] $gnt[3] $gnt[5] $gnt[4] comments: $gnt[1]\n";
	}
	else
	{
		print "$gnt[1] alleles $gnt[4] $gnt[5] not equal to reference uc($sys)\n";
	}
	
}
close(WRGNT);
close(GNT);

#open(DBSNP,$dbsnp) or die "no file $dbsnp found\n";
#while(<DBSNP>)
#{
#	chomp($_);
#	@dbsnp=split("\t",$_);
#	if(exists($gnt{$dbsnp[4]}))
#	{
#		$dbsnp{$dbsnp[4]}=$dbsnp[8];
#		print "$dbsnp[4] $dbsnp[8]\n";
#	}	
#}
